library(ggplot2)
library(patchwork)

# Study Models ------------------------------------------------------------
setwd("./Code/")
source("./main/study2_sems.R")
source("./main/study3_sems.R")
source("./main/study4_sems.R")

# Combine Studies ---------------------------------------------------------
ests_wide <- rbind(ests_wide_s2,
                   ests_wide_s3,
                   ests_wide_s4)

ests_wide$study <- factor(ests_wide$study,
                          labels = c("Study 2 (Dynata)", "Study 3 (MTurk)", "Study 4 (Bovitz)"))


# Summary Stats -----------------------------------------------------------
## Average Correlation Differences
summary(abs(ests_wide$`est.std.No Method` - 
              ests_wide$est.std.Method))
sd(abs(ests_wide$`est.std.No Method` - 
         ests_wide$est.std.Method))

# S2
mean(abs(ests_wide_s2$`est.std.No Method`-ests_wide_s2$est.std.Method))
# S3
mean(abs(ests_wide_s3$`est.std.No Method`-ests_wide_s3$est.std.Method))
# S4
mean(abs(ests_wide_s4$`est.std.No Method`-ests_wide_s4$est.std.Method))


# Proportion of corrections increasing correlation
prop.table(table(ests_wide$`est.std.No Method` <
                   ests_wide$est.std.Method))

# S2
prop.table(table(ests_wide_s2$`est.std.No Method` < ests_wide_s2$est.std.Method))
# S3
prop.table(table(ests_wide_s3$`est.std.No Method` < ests_wide_s3$est.std.Method))
# S4
prop.table(table(ests_wide_s4$`est.std.No Method` < ests_wide_s4$est.std.Method))

## FN 16
summary(abs(ests_wide$`est.std.No Method`[c(5, 7:42)] - 
              ests_wide$est.std.Method[c(5, 7:42)]))
sd(abs(ests_wide$`est.std.No Method`[c(5, 7:42)] - 
         ests_wide$est.std.Method[c(5, 7:42)]))
# Proportion of corrections increasing correlation
prop.table(table(ests_wide$`est.std.No Method`[c(5, 7:42)] <
                   ests_wide$est.std.Method[c(5, 7:42)]))


# Correlation-Specific Changes --------------------------------------------------
ests_wide$corr_diff <- abs(ests_wide$`est.std.No Method` - ests_wide$est.std.Method)

ests_wide[order(ests_wide$corr_diff, decreasing = T), 
          c("study","lhs", "rhs", "est.std.No Method", "est.std.Method", "corr_diff", "pvalue.No Method", "pvalue.Method")]


cd_hs <- subset(ests_wide, rhs == "hs" | lhs == "hs")
cd_rr <- subset(ests_wide, rhs == "rr" | lhs == "rr")
cd_pop <- subset(ests_wide, rhs == "pop" | lhs == "pop")
cd_nfc <- subset(ests_wide, rhs == "nfc" | lhs == "nfc")
cd_viol <- subset(ests_wide, rhs == "viol" | lhs == "viol")
cd_consp <- subset(ests_wide, rhs == "consp" | lhs == "consp")
cd_antdem <- subset(ests_wide, rhs == "antidem" | lhs == "antidem")

mean(cd_hs$corr_diff[which(cd_hs$study != "Study 2")]) # .07, .06 w/o S2
mean(cd_rr$corr_diff[which(cd_rr$study != "Study 2")]) # .04, .04
mean(cd_pop$corr_diff[which(cd_pop$study != "Study 2")]) # .08, .05
mean(cd_nfc$corr_diff[which(cd_nfc$study != "Study 2")]) # .04, .04
mean(cd_viol$corr_diff[which(cd_viol$study != "Study 2")]) # .02, .02
mean(cd_consp$corr_diff[which(cd_consp$study != "Study 2")]) # .06, .04
mean(cd_antdem$corr_diff[which(cd_antdem$study != "Study 2")]) # .06, .03


# Compare to Naive Corrs --------------------------------------------------
obvar_corrs <- read.csv("../Data/naive_corr_labeled.csv")

obvar_corrs$alignrev <- factor(obvar_corrs$alignrev,
                               levels = c(0, 1),
                               labels = c("PW-PW", "NW-NW"))

obvar_corrs <- obvar_corrs[which(duplicated(obvar_corrs) == FALSE), ]

obvar_corrs$study <- paste("Study", obvar_corrs$study)

ests_wide$pair <- paste(ests_wide$lhs, ests_wide$rhs, sep = "-")
ests_wide$study2 <- NA
ests_wide$study2[which(ests_wide$study == "Study 2 (Dynata)")] <- "Study 2"
ests_wide$study2[which(ests_wide$study == "Study 3 (MTurk)")] <- "Study 3"
ests_wide$study2[which(ests_wide$study == "Study 4 (Bovitz)")] <- "Study 4"

corr_compare <- merge(ests_wide, obvar_corrs, by.x = c("study2", "pair"), by.y = c("study", "pair"))

corr_compare$corr_diff <- corr_compare$est.std.Method-corr_compare$cor1


# Descriptives ------------------------------------------------------------
# Descriptives: All Scales
# correlations that increase
prop.table(table(corr_compare$corr_diff > 0))

# mean, range, SD
summary(abs(corr_compare$corr_diff))
sd(abs(corr_compare$corr_diff))


# Descriptives: By Alignment
# correlations that increase
prop.table(table(corr_compare$corr_diff[which(corr_compare$alignrev == "PW-PW")] > 0))
prop.table(table(corr_compare$corr_diff[which(corr_compare$alignrev == "NW-NW")] > 0))

# mean, range, SD
summary(abs(corr_compare$corr_diff[which(corr_compare$alignrev == "PW-PW")]))
sd(abs(corr_compare$corr_diff[which(corr_compare$alignrev == "PW-PW")]))

summary(abs(corr_compare$corr_diff[which(corr_compare$alignrev == "NW-NW")]))
sd(abs(corr_compare$corr_diff[which(corr_compare$alignrev == "NW-NW")]))


# Export to Stata for Plot ------------------------------------------------
# write.dta(corr_compare, file = "../data/mainplot.dta")


